This project set out to measure the direct and indirect effects of spatiotemporal variation in salinity on intertidal community structure in the Strait of Georgia, British Columbia. We asked the following questions:
We predicted that observed differences in community structure between consistently high salinity environments and periodically low salinity environments would be driven by trophic effects as well as physiological stress. We hypothesized that limpets and other marine gastropods are disproportionately affected by the annual decrease in surface salinity experienced in the coastal areas near the mouth of the Fraser River. We further predicted that the abundance of gastropods in periodically low salinity environments would be lower than that of consistently high salinity areas around the Southern Gulf Islands. This reduction in gastropod grazers would lead to a reduction in grazing pressure on palatable green algae, biofilms and algal spores, resulting in greater algal cover in low salinity environments. Increased algal abundance has been shown to negatively impact barnacle species (Farrell 1991), and we therefore predicted lower barnacle abundance in low salinity sites. To test these predictions, we combined laboratory salinity tolerance trials on limpets and green algae with transect surveys and limpet exclusion/inclusion experiments at three sites within each of the two salinity regions.
We conducted field studies at three sites within each of two regions with contrasting salinity regimes: West Vancouver, which experiences reduced salinities during the summer, and the Southern Gulf Islands, which experience consistently high salinities year-round. West Vancouver is located approximately 30km north of the Fraser River outflow, and the Southern Gulf Islands are located approximately 30 km southwest (Fig. 1). Gulf island sites are located on the southwest side of the island chain, and are not exposed to the Fraser River plume. The two areas are similar in terms of climate and topography. Sea surface temperature in the two regions is comparable, ranging from 5.0 to 18.5 in West Vancouver and 6.0 to 18.5 in the Gulf Islands (Fisheries and Oceans Canada 2009). The tidal range is greater in West Vancouver, with extreme high tides reaching 4.7 m above Canadian chart datum (approximated as the lowest astronomical tide), compared to 3.4 m above chart datum in the Gulf Islands. All sites were composed of granite rock except HS1, which was sandstone. The slope of the rock face ranged from 1 to 37º and aspect ranged from 40 to 320º east of magnetic north.
Figure 1: Map of the study region. Low salinity sites are located in West Vancouver and high salinity sites are located in the Gulf Islands.
Transect Surveys
Surveys were conducted once per month during low tide from May to August 2011, at each of the six study sites. Because the tidal range differs between the two areas, surveys were conducted at the vertical height corresponding to approximately 30% submersion time. This occurs at 2.1 m in the Southern Gulf Islands and 3.0 m in West Vancouver. Ten meters of transect tape were laid across the selected area and eight randomly selected points were surveyed using a 25x25 cm quadrat. Mobile invertebrates were counted and sessile invertebrates and algae were measured by percent cover.
Field Exclusion Experiment
Seven sub-sites within each of the six study sites (HS1, HS2, HS3, LS1, LS2 and LS3) were manually cleared of organisms. Within each sub-site, a limpet inclusion, exclusion and control plot were set. Inclusions and exclusions were formed by securing two copper fences, 2.5 cm high and 28.5 cm in diameter, to the rock face using Quickcrete® quick drying cement. Copper enclosures/exclosures of this type have been shown to be effective barriers to limpets (Harley 2002) and partial barriers to periwinkles (Harley 2006). Four Lottia pelta, 20±5 mm in diameter, were collected from the nearby shoreline within each site and placed into rings designated as limpet inclusion treatments. This density of 0.63 limpets per 100 cm², approximately corresponds to the average density of limpets in low salinity sites as determined by previous survey data (0.55 per 100 cm² in low salinity sites and 3.65 per 100 cm² in high salinity sites). L. pelta were determined to be the largest grazer present in the mid-intertidal in both salinity regions, and body size of grazers has been shown to positively correlate with grazing pressure (Geller 1991). Any other grazers found inside rings were removed. Limpet treatment densities were maintained every two weeks by adding or removing limpets as necessary. One circular plot within each area, also 28.5 cm in diameter, was marked with steel bolts and served as a control. The position of control and treatments was randomized within each subsite. Copper controls were not included in this study, as previous work has shown that partial copper treatments lead to partial effects which are difficult to interpret (Johnson 1992).
Sampling occurred once per month during low tide from May to August. A 10x10 cm quadrat was used to count mobile invertebrates and barnacles and estimate percent cover of algae and mussels within each treatment. Salinity samples were taken at each sampling event and measured using a refractometer (S/Mill-E, Atago Inc.).
All analyses were performed in R version 4.1.2.
Spatial and Temporal Variation in Salinity
* Salinity data were analyzed with ?
Multivariate Analyses
Univariate Analyses
Transect Surveys
Species abundance over time in low and high salinity sites were modelled with generalised linear models using the glmmTMB package version 1.1.2.3. Each species was modelled with month and salinity region as fixed factors, and a random effect of site nested within salinity region. For mobile invertebrates that were counted as total individuals (compared to percent cover), models were fit with with a negative binomial and poisson error distribution, with and without zero inflation. Sessile invertebrates and algal species were counted as percent cover of quadrats, and therefore we also fit models with a beta and a tweedie error distribution to these species. Models were then compared with the Akaike information criterion (AIC), and whichever had the lowest AIC value was chosen. Model diagnostics (residuals plots and Durbin-Watson tests) were run using the DHARMa package, version 0.4.4 (Hartig 2020). If variance across groups was unequal, we included a dispersion formula for either month, salinity region or both. Limpets, littorines, Chthamalus dalli, Red algae, Green algae and Brown algae were all fitted with a negative binomial distribution. Balanus glandula was fitted with a beta error distribution. To deal with the excessive number of zeros in the limpet count data, we also modelled these data with a zero-inflated negative binomial model. The significance of fixed effects was then tested using Wald $X^{2} tests in the car package version 3.0-11. Differences among the levels of the effects were tested using post-hoc pairwise comparisons, based on estimated marginal means using a Tukey adjustment with the emmeans package, version 1.6.3.
Field Exclusion Experiment
To test whether our exclusion plots were effective at keeping grazers out, we pooled all grazer species (limpets and littorines) and fit a generalized linear model with a negative binomial error distribution and a dispersion formula for salinity region. Because all months prior to July were likely to impact community structure, we included data from May, June and July in this model only. We fit generalized linear models to B. glandula, C. dalli, all algae pooled together, algae split into taxonomic groups of red, brown and green, as well as split into species which are grazed and non-grazed by limpets. There was minimal amounts of red and brown algae, thus we are only including models for green algae, and all algal species pooled together. C. dalli, all algae species, and green algae were all modelled as zero-inflated data as well. We were unable to fit a model including a random intercept of site nested within region to the B. glandula data without resulting in convergence issues. To resolve this, we fit the model with site as a fixed effect instead.
Spatial and Temporal Variation in Salinity
Surface salinities showed a seasonal pattern across all sites, with highest salinities occurring in the winter and early spring months (October to April) and lowest salinities occurring in the summer and early autumn (June to September) (Fig. 2). Salinity in the Gulf Islands was consistently higher than in West Vancouver (ANOVA, P<0.0001), and this difference was greatest in the summer months (May to August).
Figure 2: Measured surface salinity (psu) from sites in the Gulf Islands (shaded) and in West Vancouver (unshaded), British Columbia. Dashed line indicates Fraser River discharge rate (10³m³/s) measured at Hope, British Columbia (Environment Canada, 2012). Surface salinity for Eagle Cove, April 7, 2011, was influenced by heavy rainfall.
Transect Surveys
A stable NMDS plot with 3 dimensions and a stress of 0.15 was reached after 30 tries. PERMANOVA showed a significant effect of salinity region (F = 81.176, p = 0.005; Table 1), month (F = 4.481, p = 0.005; Table 1), and the interaction between the two (F = 1.697, p = 0.005; Table 1) on community composition.
Plots from low and high salinity sites closely cluster together. In the low salinity sites, surveys from July and August begin to form slightly less variable clusters, which may be explained by the seasonal change in salinity at these sites.
Figure 3: nMDS plot of community data collected from transect surveys from 3 high and 3 low salinity sites, collected summer 2011
| treatments | df | SS | MS | pseudo.F | R2 | p.value |
|---|---|---|---|---|---|---|
| salinity region | 1 | 14.76 | 14.76 | 81.18 | 0.29 | 0.005 |
| month | 3 | 2.44 | 0.81 | 4.48 | 0.05 | 0.005 |
| salinity region x month | 3 | 0.93 | 0.31 | 1.7 | 0.02 | 0.005 |
| residuals | 184 | 33.46 | 0.18 | 0.65 | ||
| total | 191 | 51.60 | 1.00 |
A PERMDISP test revealed that dispersion among salinity groups was equal (F = 0.23, P > 0.05; table 2), but that May and June have less variance compared to September (F = 6.515, p = 0.001; table 3). As groups with different dispersions may result in misleadingly low p values, PERMANOVA results must be interpreted cautiously.
Figure 4: dispersion levels of centroids among salinity levels and during each month of sampling
| Df | SS | MS | F | N.Perm | Pr..F. |
|---|---|---|---|---|---|
| 1 | 0.0041347 | 0.0041347 | 0.229678703154221 | 999 | 0.662 |
| 190 | 3.4204103 | 0.0180022 |
| Df | SS | MS | F | N.Perm | Pr..F. |
|---|---|---|---|---|---|
| 3 | 0.1408684 | 0.0469561 | 6.51498911371278 | 999 | 0.002 |
| 188 | 1.3549909 | 0.0072074 |
Simper analysis
Which species contribute at least 70% of the differences between groups?
| species | contribution | p_value | group |
|---|---|---|---|
| mytilus_pt | 0.1590332 | 0.115 | salinity |
| balanus_pt | 0.0862986 | 0.495 | salinity |
| fucus_pt | 0.0848790 | 0.115 | salinity |
| chthamalus_pt | 0.0811211 | 0.115 | salinity |
| mastocarpus_crust_pt | 0.0434617 | 0.585 | salinity |
| species | contribution | group | p_value |
|---|---|---|---|
| mytilus_pt | 0.1329027 | July_September | 0.025 |
| balanus_pt | 0.1067202 | July_September | 0.005 |
| fucus_pt | 0.0835693 | July_September | 0.035 |
| chthamalus_pt | 0.0544253 | July_September | 0.915 |
| ulva_pt | 0.0464845 | July_September | 0.005 |
| mytilus_pt | 0.1062865 | June_July | 0.990 |
| balanus_pt | 0.0907696 | June_July | 0.145 |
| fucus_pt | 0.0827536 | June_July | 0.070 |
| chthamalus_pt | 0.0687799 | June_July | 0.105 |
| barnacle_recruits_pt | 0.0471294 | June_July | 0.105 |
| mytilus_pt | 0.1283783 | June_September | 0.085 |
| balanus_pt | 0.0867641 | June_September | 0.340 |
| fucus_pt | 0.0797205 | June_September | 0.220 |
| chthamalus_pt | 0.0698817 | June_September | 0.055 |
| barnacle_recruits_pt | 0.0471294 | June_September | 0.105 |
| mytilus_pt | 0.1127107 | May_July | 0.845 |
| balanus_pt | 0.0945008 | May_July | 0.035 |
| fucus_pt | 0.0774119 | May_July | 0.455 |
| chthamalus_pt | 0.0531743 | May_July | 0.905 |
| pelta_no | 0.0453223 | May_July | 0.015 |
| mytilus_pt | 0.1068071 | May_June | 0.970 |
| barnacle_recruits_pt | 0.0677211 | May_June | 0.005 |
| chthamalus_pt | 0.0671705 | May_June | 0.080 |
| fucus_pt | 0.0658586 | May_June | 0.995 |
| balanus_pt | 0.0550453 | May_June | 1.000 |
| mytilus_pt | 0.1335991 | May_September | 0.005 |
| balanus_pt | 0.0853340 | May_September | 0.530 |
| fucus_pt | 0.0741613 | May_September | 0.800 |
| chthamalus_pt | 0.0522939 | May_September | 0.895 |
| ulva_pt | 0.0454040 | May_September | 0.005 |
A critique of simper analysis is that generally the species with the most variance also have the highest contribution. Even if you make groups that are copies of each other, the method will single out species with high contribution, but these are not contributions to non-existing between-group differences but to random noise variation in species abundances. Here I’m plotting species variance across all plots against their calculated contributions (Figure 5).
Figure 5: variance vs contribution to dissimilarity matrix of each species
Community Trajectories
The three low salinity sites experience dips in salinity seasonally as the Fraser River discharge increases. If salinity is a major determinant in the structure of intertidal communities, either directly through physiological tolerance of species or indirectly through changes in species interactions, we would expect larger changes in the trajectories of the three sites located closer to the FRE. In order to visualize this, we ordinated each sites trajectories in nMDS space (Figure 6). The three high salinity sites are all clustered close together and only diverge slightly throughout summer, whereas the three low salinity sites vary considerably at each time point.
Figure 6: Trajectories of intertidal communities at low and high salinity sites during the summer of 2011. Species abundances were averaged across 8 quadrats along a transect at each site and time point, and a matrix of Bray-Curtis dissimilarities was consructed.
Field Exclusion Experiment
Species abundance: A stable NMDS plot with 3 dimensions and a stress of 0.11 was reached after 20 tries. PERMANOVA analysis revealed significant differences in the community composition of salinity region, grazer treatment and the interaction between the two.
We hypothesized that the exclusion plots in the high salinity sites would resemble the community structure in both the control and exclusion plots of the low salinity sites, due to the lack of limpet grazing. Both treatment (pseudo-F = , p value = 0.005; Table 6), as well as the interaction between salinity region and treatment (pseudo-F, p value = 0.025; Table 6) were significant.
Figure 7: nMDS plot of community data from field exclusion experiment in high and low salinity sites from the summer of 2011.
| treatments | df | SS | MS | pseudo.F | R2 | p.value |
|---|---|---|---|---|---|---|
| salinity region | 1 | 0.97 | 0.97 | 8.2 | 0.09 | 0.375 |
| treatment | 1 | 0.40 | 0.4 | 3.41 | 0.04 | 0.005 |
| salinity region x treatment | 1 | 0.18 | 0.18 | 1.52 | 0.02 | 0.025 |
| residuals | 80 | 9.43 | 0.12 | 0.86 | ||
| total | 83 | 10.97 | 1.00 |
A PERMDISP test revealed that dispersion among salinity groups was unequal (F = 17.3, P = 0.001; table 6), as well as dispersion among treatment groups (F = 7.5, p = 0.007; table 7). As groups with different dispersions may result in misleadingly low p values, PERMANOVA results must be interpreted cautiously.
Figure 8: dispersion levels of centroids among salinity levels and treatment levels
| Df | SS | MS | F | N.Perm | Pr..F. |
|---|---|---|---|---|---|
| 1 | 0.4511003 | 0.4511003 | 17.3184440175554 | 999 | 0.001 |
| 82 | 2.1358863 | 0.0260474 |
| Df | SS | MS | F | N.Perm | Pr..F. |
|---|---|---|---|---|---|
| 1 | 0.1553742 | 0.1553742 | 7.46326674014969 | 999 | 0.008 |
| 82 | 1.7071186 | 0.0208185 |
Simper analysis
How much does each species contribute to the overall dissimilarity matrix among salinity regions (Table 8), and treatments (Table 9).
| species | contribution | p_value | group |
|---|---|---|---|
| ulva_pt | 0.2125875 | 0.115 | salinity |
| balanus_no | 0.1133128 | 0.390 | salinity |
| chthamalus_no | 0.0857018 | 0.305 | salinity |
| masto_crust_pt | 0.0630440 | 0.685 | salinity |
| fucus_pt | 0.0556836 | 0.115 | salinity |
| species | contribution | p_value | group |
|---|---|---|---|
| ulva_pt | 0.1927748 | 0.005 | treatment |
| balanus_no | 0.1076078 | 0.010 | treatment |
| chthamalus_no | 0.0872806 | 0.005 | treatment |
| masto_crust_pt | 0.0622233 | 0.785 | treatment |
| fucus_pt | 0.0546080 | 0.095 | treatment |
Transect Surveys
Invertebrates
All limpet species were pooled and include primarily Lottia paradigitalis, but also Lottia digitalis, Lottia pelta, Lottia scutum, and Lottia persona. Limpet abundance is consistently higher in the high salinity region, and saw a decreased over the course of the summer in the low salinity region only. Littorines were also pooled together and consisted of Littorina sitkana and Littorina scutulata. Littorine abundance was only higher in the high salinity region during the month of July, before decreasing slightly; Littorines also experienced a small dip in abundance during July in the low salinity region. Abundance of Balanus glandula increased over the summer in both regions, and Chthamalus dalli was consistently higher in the high salinity sites (Figure 9).
Figure 9: Invertebrate species surveyed in 3 high and 3 low salinity sites during summer 2011
| source | statistic | df | p.value |
|---|---|---|---|
| limpets | |||
| month | 22.2133596 | 3 | 0.0000589 |
| salinity_regime | 13.2517681 | 1 | 0.0002723 |
| month:salinity_regime | 27.5102771 | 3 | 0.0000046 |
| littorines | |||
| month | 26.7784889 | 3 | 0.0000066 |
| salinity_regime | 0.0009192 | 1 | 0.9758133 |
| month:salinity_regime | 27.9464824 | 3 | 0.0000037 |
| balanus | |||
| month | 41.1194879 | 3 | 0.0000000 |
| salinity_regime | 1.2372062 | 1 | 0.2660102 |
| chthamalus | |||
| month | 4.9897982 | 3 | 0.1725457 |
| salinity_regime | 66.5440393 | 1 | 0.0000000 |
Algae by taxonomic group
Red Algae consists of primarily Pyropia, and Mastocarpus (foliose and crust form), but also Hildenbrandia, Endocladia, Microcladia, and Polysiphonia. Green algae consists of primarily Ulva spp, but also Urospora, and Acrosiphonia. Brown algae consists primarly of Fucus distichus, but also Melanosiphon spp, and Leathesia spp. There is weak evidence of a decrease in red algae over the summer in both regions (F= 6.26, p = 0.09; table 11), but was consistently higher in the high salinity region (F = 13.25, p < 0.001; table 11). Green and brown algae increased in both regions over the summer (Figure 10).
Figure 10: Algae divided into taxonomic groups, surveyed at 3 high and 3 low salinity sites during the summer of 2011.
| source | statistic | df | p.value |
|---|---|---|---|
| brown algae | |||
| month | 15.8875372 | 3 | 0.0011958 |
| salinity_regime | 1.3204216 | 1 | 0.2505164 |
| green algae | |||
| month | 37.8466010 | 3 | 0.0000000 |
| salinity_regime | 0.4019598 | 1 | 0.5260789 |
| red algae | |||
| month | 6.2589556 | 3 | 0.0996691 |
| salinity_regime | 13.2491129 | 1 | 0.0002727 |
Field Exclusion Experiment
Grazersfigure 11: Grazer abundance in grazer exclusion vs. control plots in high and low salinity sites.
| source | statistic | df | p.value |
|---|---|---|---|
| treatment | 0.0612935 | 1 | 0.8044631 |
| region | 10.0465690 | 1 | 0.0015263 |
| treatment:region | 7.7226311 | 1 | 0.0054533 |
Barnacles
There were site level differences in B. glandula (F = 65.31, p < 0.001; Fig.12), and weak evidence of a treatment level effect (F = 2.64, p = 0.10; Fig.12) . C. dalli decreased in the grazer exclusion plots in the high salinity sites only.Figure 12: Barnacle abundance in grazer exclusion vs. control plots in high and low salinity sites.
| source | statistic | df | p.value |
|---|---|---|---|
| treatment | 2.640882 | 1 | 0.1041457 |
| site | 65.318524 | 5 | 0.0000000 |
| treatment:site | 8.834319 | 5 | 0.1158579 |
| source | statistic | df | p.value |
|---|---|---|---|
| treatment | 0.0477698 | 1 | 0.8269905 |
| region | 2.7900916 | 1 | 0.0948488 |
| treatment:region | 17.1414786 | 1 | 0.0000347 |
All algae
Algae percent cover increased in the grazer exclusion plots in the high salinity sites (figure 13). This includes reds, greens, browns and diatoms.Figure 13: Algae abundance in grazer exclusion vs. control plots in high and low salinity sites.
| source | statistic | df | p.value |
|---|---|---|---|
| treatment | 1.243033 | 1 | 0.2648873 |
| region | 19.010591 | 1 | 0.0000130 |
| treatment:region | 19.127575 | 1 | 0.0000122 |
Figure 14: Green algae abundance in grazer exclusion vs. control plots in high and low salinity sites.
| term | statistic | df | p.value |
|---|---|---|---|
| treatment | 0.6632273 | 1 | 0.4154229 |
| region | 18.6691845 | 1 | 0.0000155 |
| treatment:region | 17.5541920 | 1 | 0.0000279 |
Because plots were cleared in May, species were unlikely to have reached carrying capacity by July, and species abundance may not accurately depict the impact of the experimental manipulation. To account for this in our analysis, we also collapsed the data to presence/absence and re-ran the ordination and permutation as described above, on a matrix of Jaccard dissimilarities.
Figure S1: nMDS plot of species richness data from field exlcusion experiment in high and low salinity sites.
A PERMANOVA on community composition of the presence or absence of species found treatment (pseudo F = 8.97, p value = 0.005; Table S1), salinity region (pseudo F = 19.8, p value = 0.005; Table S1) and the interaction between the two (pseudo F = 8.38, p value = 0.005; Table S1) were significant.
| treatments | df | SS | MS | pseudo.F | R2 | p.value |
|---|---|---|---|---|---|---|
| salinity region | 1 | 2.28 | 2.28 | 19.8 | 0.17 | 0.005 |
| treatment | 1 | 1.03 | 1.03 | 8.97 | 0.08 | 0.005 |
| salinity region x treatment | 1 | 0.97 | 0.97 | 8.38 | 0.07 | 0.005 |
| residuals | 80 | 9.22 | 0.12 | 0.68 | ||
| total | 83 | 13.51 | 1.00 |
Algae by functional group
We also grouped algae by functional groups. ‘Grazed’ algae are filaments and sheets that are more likely to be consumed by invertebrates at a higher frequency, and consist of Ulva, Pyropia, Polysiphonia, Diatoms, Urospora, Acrosiphonia. ‘Non-grazed’ algae are turfs and crusts that are consumed at a lower frequency, and consist of Fucus, Mastocarpus, Mastocarpus crust, Hildenbrandia, Endocladia, Melanosiphon, Leathesia. Grazed algae increased over the summer in the low salinity region only, and non-grazed algae increased over the summer in both regions (Figure S2).
Figure S2: Algae divided into functional groups, surveyed at 3 high and 3 low salinity sites during the summer of 2011.
| source | statistic | df | p.value |
|---|---|---|---|
| grazed algae | |||
| month | 18.0978993 | 3 | 0.0004199 |
| salinity_regime | 2.5636368 | 1 | 0.1093468 |
| month:salinity_regime | 9.8419923 | 3 | 0.0199581 |
| non grazed algae | |||
| month | 17.2886876 | 3 | 0.0006164 |
| salinity_regime | 0.7150894 | 1 | 0.3977594 |
| month:salinity_regime | 3.4599121 | 3 | 0.3259987 |
Figure S3: Abundance of algae divided by taxonomic groups in grazer exclusion vs. control plots in high and low salinity sites.
| source | statistic | df | p.value |
|---|---|---|---|
| brown algae | |||
| treatment | 1.7359146 | 1 | 0.1876573 |
| region | 2.3617753 | 1 | 0.1243406 |
| treatment:region | 0.0000012 | 1 | 0.9991184 |
| green algae | |||
| treatment | 0.6632273 | 1 | 0.4154229 |
| region | 18.6691845 | 1 | 0.0000155 |
| treatment:region | 17.5541920 | 1 | 0.0000279 |
| red algae | |||
| treatment | 0.4335256 | 1 | 0.5102640 |
| region | 0.2108588 | 1 | 0.6460951 |
| treatment:region | 0.1020306 | 1 | 0.7494063 |
Figure S4: Abundance of algae divided by functional group in grazer exclusion vs. control plots in high and low salinity sites.
| source | statistic | df | p.value |
|---|---|---|---|
| grazed algae | |||
| treatment | 0.5250043 | 1 | 0.4687148 |
| region | 1.5083479 | 1 | 0.2193914 |
| treatment:region | 2.5625582 | 1 | 0.1094214 |
| non grazed algae | |||
| treatment | 1.4517840 | 1 | 0.2282419 |
| region | 0.8582985 | 1 | 0.3542153 |
| treatment:region | 0.3710508 | 1 | 0.5424320 |